Multi-cluster and Traveling Chimera States in Phase-Coupled Oscillators

Hsien-Ching (Kao^1), Jianbo (Xie^2) & Edgar (Knobloch^3)

(( )^1)WTG Group, Wolfram Research, Inc., Champaign, IL
(( )^(2,3))Department of Physics, University of California, Berkeley, CA

hkao@wolfram.com

October 24, 2014

Coupled Oscillators Systems

Arrays or networks of limit-cycle oscillators that couple locally or globally.

Physics: laser arrays, coupled Josephson junctions

Biology: firefly population, neuron systems

Engineering: power grids, mechanical oscillators arrays

These systems are frequently modeled by

Complex Ginzburg-Landau equations

(Subscript[A, t])= μA+(1+iα)(Subscript[A, xx])-(1+iβ)|A((|)^2)A+f(A,(SuperStar[A]))

Reaction diffusion systems

(Subscript[u, t])= M (Subscript[u, xx])+N(u), u∈(R^n)

Phase Coupled Oscillators

Limit cycle (Fitzhugh-Nagumo model)

v'=v-((v^3)/3)-w+0.5,  13 w'=v+0.7-0.8w

Animate[
Row[{
ParametricPlot[#,{t,0,T},PlotRange->{{-2.2,2.2},{-0.5,1.7}},PlotPoints->100,PlotStyle->Black,Frame->True, FrameLabel->{Style["v",20],Style["w",20]},ImageSize->{250,250}],Spacer[25],
Plot[#,{t,0,T},PlotRange->{{0,200},{-2.0,2.0}},ImageSize->{250,250},PlotStyle->{Black,Red},Frame->True, FrameLabel->{Style["t",20],None},PlotLegends->{"v","w"}]
}]&[{v[t],w[t]}/.sol1],{{T,10,"t"},1,200},AnimationRunning->False,Initialization:>(sol1=NDSolve[{v'[t]==v[t]-(((v[t])^3)/3)-w[t]+0.5,13 w'[t]==v[t]+0.7-0.8 w[t],v[0]==1,w[0]==0},{v,w},{t,0,200}]),TrackedSymbols:>T]

Phase reduction (Malkin’s theorem)

((d (Subscript[X, i]))/(d t))=(Subscript[F, i])((Subscript[X, i]))+ϵ (Subscript[p, i])((Subscript[X, 1]),(Subscript[X, 2,])...,(Subscript[X, i]),...)

The system can be approximated as

((d (Subscript[θ, i]))/(d t))=(Subscript[ω, i])+ϵ (Subscript[Γ, i])((Subscript[θ, 1]),(Subscript[θ, 2]),...,(Subscript[θ, i,])...)

Chimera State

Kuramoto and Battogtokh (2002)

(Subscript[θ, t])=ω - (∫_(-π)^π)G(x-y)Sin[θ(x,t)-θ(y,t)+α]dy,

where G(x) is proportional to exp(-κ|x|).

Chimera state (partially synchronized in space)

n=256;
κ=4.0;
α=1.457;
dt=0.025;
dtfactor=2.Pi dt;
rk4coef={1.,2.,2.,1.}/6;
fourier=Fourier[#,FourierParameters->{-1, -1}]&;
ifourier=InverseFourier[#,FourierParameters->{-1, -1}]&;
xx=N[2.Pi Range[0,n-1]/n];
With[{xx1= Min[#,2Pi-#]&/@Abs[xx]},
kernel=fourier[Exp[-κ/(2Pi)xx1]];
phiinitial=Exp[-(30/((2Pi)^2))xx1^2]RandomReal[{-3.,3.},n];];
rhs=(-dtfactor)Im[Exp[I(#+α)]ifourier[kernel*fourier[Exp[-I#]]]]&;
timestepping[sol_List]:=sol+rk4coef.FoldList[rhs[sol+#1 #2] &,rhs[sol],{0.5,0.5,1.}];
kuramoto=With[{phitemp=Nest[timestepping,phiinitial,7000]},
NestList[timestepping,phitemp,1000]];
Save["kuramoto.txt",kuramoto];

Panel[Row[{,Spacer[50.],Animate[ListPlot[Thread[{xx1,Mod[kuramoto[[index]]+1.5-m index,2Pi]-Pi}],PlotStyle->Black,PlotRange->{{0,2Pi},{-Pi,Pi}},AxesLabel->{Style["x",15],Style["θ",15]},ImageSize->250],
{{index,1,Style["",15]},1,Length[kuramoto],1},
AnimationRunning->False,
TrackedSymbols:>index,
Initialization:>(
Clear[kuramoto];
n=256;
xx1=N[2.Pi Range[0,n-1]/n];
Get["kuramoto.txt"];
m=Last@LinearModelFit[kuramoto[[All,128]],x,x]["BestFitParameters"];
)]},Alignment->{Center,Center}]]

Modeling Equation

(Subscript[θ, t])=ω - (∫_(-π)^π)G(x-y)Sin[θ(x,t)-θ(y,t)+α]dy, α∈[0,π/2],

with G(x) equal to (Subscript[G, n]^((1)))(x)=Cos(n x) and (Subscript[G, n]^((2)))(x)=Cos(n x)+Cos((n+1)x). ω can be set to 0 without losing of generality.

Effective equation

(Subscript[z, t])=(1/2)[(e^(-i α))Z(x,t)-(z^2)(e^(i α))(SuperStar[Z])(x,t)], z(x,t)≡exp[i θ(x,t)],

where Z=K[z] and the linear operator K is defined as

K[u](x,t)≡(∫_(-π)^π)G(x-y)u(y,t)dy

Rotating wave solution

z(x,t)=(e^(-i Ω t))(OverTilde[z])(x)

Using this ansatz with the effective equation one can derive the self-consistency equation

(OverTilde[Z])(x)≡R(x)(e^(i Θ(x)))=(∫_(-π)^π)G(x-y)(e^(i β))((Ω-μ(y))/((SuperStar[(OverTilde[Z])])(y)))dy,

where μ(x)≡({{\[Piecewise], ({{((((Ω^2)-((\[LeftBracketingBar](OverTilde[Z])(x)\[RightBracketingBar])^2)))^(1/2)), ,Ω≥\[LeftBracketingBar](OverTilde[Z])(x)\[RightBracketingBar] }, {((i(((\[LeftBracketingBar](OverTilde[Z])(x)\[RightBracketingBar])^2)-(Ω^2)))^(1/2)), ,Ω<\[LeftBracketingBar](OverTilde[Z])(x)\[RightBracketingBar]}})}}), and β≡(π/2)-α.

(OverTilde[Z])≡K[(OverTilde[z])] is referred as the local order parameter.

Splay State

θ(x,t)=q x+Ω t, q∈Z,

where Ω=-(∫_(-π)^π)G(y)Sin(q y+α)dy. When α<π/2, splay states with |q|=n (G=(Subscript[G, n]^((1)))) and |q|=n, n+1 (G=(Subscript[G, n]^((2)))) are stable, respectively.

xx=2.Pi Range[512]/512;
Animate[Row[{ListPlot[Thread[{xx,splay1[[n]]}],PlotStyle->Black,PlotRange->{{0,2Pi},{0,2Pi}},ImageSize->250,AxesLabel->{Style["x",15],Style["θ",15]}],Spacer[50],ListPlot[Thread[{xx,splay2[[n]]}],PlotStyle->Black,PlotRange->{{0,2Pi},{0,2Pi}},ImageSize->250,AxesLabel->{Style["x",15],Style["θ",15]}]}],
{{n,1,Style["",15]},1,Length[splay1],1},
AnimationRunning->False,TrackedSymbols:>n,
Initialization:>(
splay1=PhaseOscillators1[2.Pi Range[512]/512,0.1,Range[0,10,0.1],1,0.01];splay2=PhaseOscillators2[4.Pi Range[512]/512,0.1,Range[0,10,0.1],2,0.01];
)]

Splay states with (Left) (Subscript[G, 1]^((1)))(x) and (Right) (Subscript[G, 2]^((2)))(x) coupling, β=0.1 and N=512

Chimera State G=(Subscript[G, n]^((1)))

Panel@Row@{Animate[ListPlot[Thread[{xx,Mod[chimera1sol[[n]]-1-m1 n,2Pi]-Pi}],PlotStyle->Black,ImageSize->300,PlotRange->{{0,2Pi},{-Pi,Pi}},AxesLabel->{Style["x",15],Style["θ",15]}],
{{n,1,Style["",15]},1,Length[chimera1sol],1},
AnimationRunning->False,TrackedSymbols:>n,
Initialization:>(
Get["chimera1.txt"];
chimera1sol=PhaseOscillators1[chimera1,0.1,Range[0,10,0.1],1,0.01];
m1=Last@LinearModelFit[chimera1sol[[;;20,128]],x,x]["BestFitParameters"];
Clear[chimera1];)],Spacer[20],With[{Z=(2.Pi/512)(Cos[xx](Cos[xx].#)+Sin[xx](Sin[xx].#))&[Exp[I First[chimera1sol]]]},GraphicsColumn[{ListPlot[Abs[Z],ImageSize->250,AspectRatio->0.5,PlotStyle->Black,AxesLabel->{Style["x",15],Style["R(x)",15]}],ListPlot[Arg[Z],ImageSize->250,AspectRatio->0.5,PlotStyle->Red,AxesLabel->{Style["x",15],Style["Θ(x)",15]}]}]]}

Chimera states with (Subscript[G, 1]^((1))) coupling, β=0.1 and N=512

(OverTilde[Z])(x) is of the form (Subscript[R, 0])cos(n x)

(Subscript[R, 0]^2)=(e^(i β ))((∫_0)^(2π)){Ω-(([(Ω^2)-(Subscript[R, 0]^2) (cos^2)(ny)])^(1/2))}dy

The integral is identical for all integer (and half integer) n.

Integral-Algebraic Equation

(R^2)=(e^(iβ))((∫_0)^(2π))Ω-(Sqrt[(Ω^2)-(R^2)(cos^2)(y)])dy

Clear[β,Ω,R,eqs,int,y];
eqs[β_]={R^2==Re[E^(I β)int[Ω,R]],0==Im[E^(I β)int[Ω,R]]};
int[Ω_?NumericQ,R_?NumericQ]:=4Total[NIntegrate[Ω-(Sqrt[Ω^2-R^2 Cos[y]^2]),#]&/@{{y,0,ArcCos[Ω/R]},{y,ArcCos[Ω/R],Pi/2}}]

When β=0, the equation admits a solution R=Ω=2π-4.

Check the solution

{sol=FindRoot[eqs[0.],{{R,2Pi-4},{Ω,2Pi-4}}],N[2Pi-4]}

Determine the dependence of R and Ω on β by pseudo-arclength continuation

f(x,α)=0, x∈(R^n), α∈R

Parametrized x and α with respect to arclength, i.e., x=x(s), α=α(s)

f(x(s),α(s))=0, x'(((s))^2)+α'(((s))^2)=1

In pseudo-arclength continuation

f(x(s),α(s))=0, x'(s)(x(s)-(Subscript[x, 0]))+α'(s)(α(s)-(Subscript[α, 0]))=Δs

sol0=Append[sol[[All,2]],0.];
sol1=Append[FindRoot[eqs[Pi/200.],{{R,2Pi-4},{Ω,2Pi-4}}][[All,2]],Pi/200.];
ds=0.1;
data=NestWhileList[(vec=Normalize[#-sol0];sol0=#;
FindRoot[{Sequence@@eqs[β],vec.({R,Ω,β}-#)==ds},Thread[{{R,Ω,β},#+vec ds}]])[[All,2]]&,sol1,(Last[#]<Pi/2.)&];

With[{R=Interpolation[data[[All,{3,1}]]],Ω=Interpolation[data[[All,{3,2}]]]},
Plot[{R[β],Ω[β]},{β,0,Pi/2},PlotStyle->{Black,Red},PlotLegends->{"R","Ω"},AxesLabel->{Style["β",20],Automatic}]//Quiet]

Stability Analysis

Linearizing the effective equation about (OverTilde[z])(x) with v(x,t)

(Subscript[v, t])=i μ v+(1/2)[(e^(-i α))V(x,t)-((OverTilde[z])^2)(e^(i α))(SuperStar[V])(x,t)], V=K[v]

Solving the equation by v(x,t)=(e^(λ t))(Subscript[v, 1])(x)+(e^((SuperStar[λ])t))(Subscript[v, 2])(x) gives

λ(({{(Subscript[v, 1])}, {(Subscript[v, 2])}}))=(1/2)(({{2 i μ+(e^(-i α))K, -(e^(i α))((OverTilde[z])^2)K}, {-(e^(-i α))((OverTilde[z])^( *2))K, -2 i (SuperStar[μ])+(e^(i α))K}}))(({{(Subscript[v, 1])}, {(Subscript[v, 2])}}))

Essential (Continuous) Spectrum: (Subscript[{i μ(x),-i (SuperStar[μ])(x)}, x∈[-π,π]])

Point spectrum

data=Partition[ReadList["eig1.txt",Number],3];
initialize=Function[{β,a,Ω},
α=π/2-β;xs=ArcCos[Ω/a];
μ=Sqrt[(Ω^2)-(((a Cos[#]))^2)]&;
z2=(((Exp[I β]a Cos[#]/(Ω+μ[#])))^2)&;
M0=DiagonalMatrix[2Exp[I α #]&/@N[{1,1,-1,-1}]];
nint=(NIntegrate[#[x],{x,0,xs},Method-> "DoubleExponential"]+NIntegrate[#[x],{x,xs,π/2},Method-> "DoubleExponential"])&;
f11=((-1)/(#2-I μ[#1]))&;f12=((z2[#1])/((#2-I μ[#1])))&;
f21[x_,λ_]:=Conjugate[z2[x]/(Conjugate[λ]-I μ[x])];
f22[x_,λ_]:=Conjugate[-1/(Conjugate[λ]-I μ[x])];
];
detM[λ_]:=
Module[{f110,f11c,f120,f12c,f210,f21c,f220,f22c},
{f110,f120,f210,f220}=nint[#]&/@{f11[#,λ]&,f12[#,λ]&,f21[#,λ]&,f22[#,λ]&};
{f11c,f12c,f21c,f22c}=
nint[#]&/@{(f11[#,λ]Cos[2#])&,(f12[#,λ]Cos[2#])&,(f21[#,λ]Cos[2#])&,(f22[#,λ]Cos[2#])&};
Det[M0+
2.0{{f110 ,f11c,f120,f12c},
{f11c,f110,f12c,f120},
{f210,f21c,f220,f22c},
{f21c,f210,f22c,f220}}]
]
ddetM=Function[{λ,h},(0.25/h)ParallelMap[detM,λ+h{1,-1,I,-I}].{1,-1,-I,I}];
newtonM[xi_,h_,tol_]:=FixedPoint[(#-((detM[#])/(ddetM[#,h])))&,xi,SameTest->(Abs[(#1-#2)/#2]<tol&)]
contourdetM=Function[{x1,x2,y1,y2,nn},
Module[{z,data},
z={x1+I Range[y1,y2,(y2-y1)/nn],I y2+Range[x1,x2,(x2-x1)/nn],
x2+I Range[y2,y1,(y1-y2)/(nn)],I y1+Range[x2,x1,(x1-x2)/nn]};
data=ParallelMap[detM,z,{2}];
-Total[Map[(#[[2;;]]/#[[1;;-2]]-#[[1;;-2]]/#[[2;;]])&,data],2]/(4I Pi)
]];

Do[If[n<68,AppendTo[data[[n]],0.],
{β,a,Ω}=data[[n]];initialize[β,a,Ω];
AppendTo[data[[n]],newtonM[1.,0.001,N[1*^-4]]]],{n,Length@data}];
Export["eigpdata1.txt",data]

data=Partition[ReadList["eigpdata1.txt",Number],4];
Panel@Row[{GraphicsColumn[{ListPlot[data[[All,{1,4}]],PlotStyle->Black,AxesLabel->{Style["β",15],Style["(Subscript[λ, p])",15]}],
Block[{β,R,Ω,λ},{β,R,Ω,λ}=data[[199]];
ListLinePlot[{{{0,Ω},{0,-Ω}},{{-Sqrt[R^2-Ω^2],0},{0,0}}},PlotStyle->Black,PlotRange->{{-3.5,1},{-2,2}},AxesLabel->{Style["Re[λ]",15],Style["Im[λ]",15]},Epilog->{PointSize[Large],Point[{λ,0}]}]]},ImageSize->200],Spacer[50],
vec1=Import["eigv1.h5",{"Datasets","/Dataset1"}];
xxx=2.Pi Range[1000]/1000;
GraphicsGrid[{{ListLinePlot[Thread[{xxx,Norm/@vec1[[1;;(Length@vec1/2)]]}],PlotRange-> Full,PlotStyle->Black,AxesLabel->{Style["x",15],Style["|(Subscript[v, 1])|",15]},Ticks->{Automatic,None}],ListLinePlot[Thread[{xxx,ArcTan@@@vec1[[1;;(Length@vec1/2)]]}],PlotRange-> Full,PlotStyle->Black,AxesLabel->{Style["x",15],Style["arg((Subscript[v, 1]))",15]}]},{ListLinePlot[Thread[{xxx,Norm/@vec1[[-(Length@vec1/2);;]]}],PlotRange-> Full,PlotStyle->Black,AxesLabel->{Style["x",15],Style["|(Subscript[v, 2])|",15]},Ticks->{Automatic,None}],ListLinePlot[Thread[{xxx,ArcTan@@@vec1[[-(Length@vec1/2);;]]}],PlotRange-> Full,PlotStyle->Black,AxesLabel->{Style["x",15],Style["arg((Subscript[v, 2]))",15]}]}},ImageSize->400]
}]

Left Top: dependence of (Subscript[λ, p]) on β
Left bottom: spectrum of the linearized operator when β=0.83
Right: eigenvectors at the onset of instability

Chimera State G=(Subscript[G, n]^((2)))

Panel@Row@{Animate[ListPlot[Thread[{xx,Mod[chimera2sol[[n]]-m2 n,2Pi]-Pi}],PlotStyle->Black,PlotRange->{{0,2Pi},{-Pi,Pi}},ImageSize->300,AxesLabel->{Style["x",15],Style["θ",15]}],
{{n,1,Style["",15]},1,Length[chimera2sol],1},
AnimationRunning->False,TrackedSymbols:>n,
Initialization:>(
Get["chimera2.txt"];
chimera2sol=PhaseOscillators2[chimera2,0.03,Range[0,10,0.1],1,0.01];
m2=Last@LinearModelFit[chimera2sol[[;;15,80]],x,x]["BestFitParameters"];
Clear[chimera2];)],Spacer[50],With[{Z=(2.Pi/512)(Cos[xx](Cos[xx].#)+Cos[2xx](Cos[2xx].#)+Sin[xx](Sin[xx].#)+Sin[2xx](Sin[2xx].#))&[Exp[I First[chimera2sol]]]},GraphicsColumn[{ListPlot[Abs[Z],ImageSize->250,AspectRatio->0.5,PlotStyle->Black,AxesLabel->{Style["x",15],Style["R(x)",15]}],ListPlot[Arg[Z],ImageSize->250,AspectRatio->0.5,PlotStyle->Red,AxesLabel->{Style["x",15],Style["Θ(x)",15]}]}]]}

Chimera states with (Subscript[G, 1]^((2))) coupling, β=0.03 and N=512

(OverTilde[Z])(x) is of the form (Subscript[R, 0])cos[(2n+1) x/2](e^(i x/2))

(Subscript[R, 0]^2)=(e^(i β ))((∫_0)^(2π)){Ω-(([(Ω^2)-(Subscript[R, 0]^2) (cos^2)((2n+1)y/2)])^(1/2))}dy

Point spectrum

data=First@Import["eigpdata2.h5","Data"];
Panel@Row[{GraphicsColumn[{ListLinePlot[{data[[All,{1,4}]],data[[All,{1,5}]],data[[All,{1,6}]]},PlotStyle->{Black,Red,Directive[Red,Dashed]},AxesLabel->{Style["β",15],Style["(Subscript[λ, p])",15]}],
Block[{β,R,Ω,λ,λr,λi},{β,R,Ω,λ,λr,λi}=data[[199]];
ListLinePlot[{{{0,Ω},{0,-Ω}},{{-Sqrt[R^2-Ω^2],0},{0,0}}},PlotStyle->Black,PlotRange->{{-3.5,1},{-2,2}},AxesLabel->{Style["Re[λ]",15],Style["Im[λ]",15]},Epilog->{PointSize[Large],Point[{λ,0}],Point[{λr,λi}],Point[{λr,-λi}]}]]},ImageSize->200],Spacer[50],
vec2=First@Import["eigv2.h5","Data"];
xxx=2.Pi Range[1500]/1500;
GraphicsGrid[{{ListLinePlot[Thread[{xxx,Norm/@vec2[[1;;(Length@vec2/2)]]}],PlotStyle->Black,AxesLabel->{Style["x",15],Style["|(Subscript[v, 1])|",15]},Ticks->{Automatic,None}],ListLinePlot[Thread[{xxx,ArcTan@@@vec2[[1;;(Length@vec2/2)]]}],PlotRange-> Full,PlotStyle->Black,AxesLabel->{Style["x",15],Style["arg((Subscript[v, 1]))",15]}]},{ListLinePlot[Thread[{xxx,Norm/@vec2[[-(Length@vec2/2);;]]}],PlotStyle->Black,AxesLabel->{Style["x",15],Style["|(Subscript[v, 2])|",15]},Ticks->{Automatic,None}],ListLinePlot[Thread[{xxx,ArcTan@@@vec2[[-(Length@vec2/2);;]]}],PlotRange-> Full,PlotStyle->Black,AxesLabel->{Style["x",15],Style["arg((Subscript[v, 2]))",15]}]}},ImageSize->400]
}]

Left Top: dependence of (Subscript[λ, p]) on β
Left bottom: spectrum of the linearized operator when β=0.83
Right: eigenvectors at the onset of instability

Chimera State G=(Subscript[G, n]^((2)))

Panel@Row@{Animate[ListPlot[Thread[{xx,Mod[chimera3sol[[n]]+0.5-m3 n,2Pi]-Pi}],PlotStyle->Black,ImageSize->300,PlotRange->{{0,2Pi},{-Pi,Pi}},AxesLabel->{Style["x",15],Style["θ",15]}],
{{n,1,Style["",20]},1,Length[chimera3sol],1},
AnimationRunning->False,TrackedSymbols:>n,
Initialization:>(
Get["chimera3.txt"];
chimera3sol=PhaseOscillators2[chimera3,0.03,Range[0,10,0.1],1,0.01];
m3=Last@LinearModelFit[chimera3sol[[;;20,150]],x,x]["BestFitParameters"];
Clear[chimera3];)],Spacer[50],With[{Z=(2.Pi/512)(Cos[xx](Cos[xx].#)+Cos[2xx](Cos[2xx].#)+Sin[xx](Sin[xx].#)+Sin[2xx](Sin[2xx].#))&[Exp[I First[chimera3sol]]]},GraphicsColumn[{ListPlot[Abs[Z],ImageSize->250,AspectRatio->0.5,PlotStyle->Black,AxesLabel->{Style["x",15],Style["R(x)",15]}],ListPlot[Arg[Z],ImageSize->250,AspectRatio->0.5,PlotStyle->Red,AxesLabel->{Style["x",15],Style["Θ(x)",15]}]}]]}

Chimera states with (Subscript[G, 1]^((2))) coupling, β=0.03 and N=512

(OverTilde[Z])(x) is of the form b sin(x)+d sin(2x), where b is real and d is complex

b (e^(-i β))=((∫_0)^(2π))((b (sin^2)y+d sin y sin 2y)/(Ω+μ(y)))dy
d (e^(-i β))=((∫_0)^(2π))((b sin y sin 2y+d (sin^2) y)/(Ω+μ(y)))dy    (Legended[, (LineLegend["\!\(\*SubscriptBox[\(b\), \(r\)]\)", "\!\(\*SubscriptBox[\(d\), \(r\)]\)", "\!\(\*SubscriptBox[\(d\), \(i\)]\)", "Ω"])])

data3=First@Import["data3.h5","Data"];
ListLinePlot[Evaluate[data3[[All,{1,#}]]&/@Range[2,5]],PlotStyle->{Black,Red,Directive[Red,Dashed],Directive[Black,Dashed]},AxesLabel->{Style["β",20],None},PlotLegends->{Style["(Subscript[b, r])",20],Style["(Subscript[d, r])",20],Style["(Subscript[d, i])",20],Style["Ω",20]}]

Coherent Traveling Wave

Manipulate[
Animate[GraphicsRow[{ListPlot[Thread[{xx,Mod[coherent2sol[[n]]+2+sign 0.199(n-1),2Pi]-Pi}],PlotStyle->Black,PlotRange->{{0,2Pi},{-Pi,Pi}},AxesLabel->{Style["x",15],Style["θ",15]}],ListPlot[Thread[{xx,Mod[coherent3sol[[n]]-1+sign 0.205(n-1),2Pi]-Pi}],PlotStyle->Black,PlotRange->{{0,2Pi},{-Pi,Pi}},AxesLabel->{Style["x",15],Style["θ",15]}]},ImageSize->600],
{{n,1,Style["",20]},1,Length[coherent2sol],1},
AnimationRunning->False,TrackedSymbols:>n,
Initialization:>(
Get["coherent2.txt"];
coherent2sol=PhaseOscillators2[coherent2,0.76,Range[0,30,0.1],1,0.01];
Clear[coherent2];
Get["coherent3.txt"];
coherent3sol=PhaseOscillators2[coherent3,0.7,Range[0,30,0.1],1,0.01];
Clear[coherent3];
)],{{sign,0,Style["Rotating Frame",15]},{0,1}},TrackedSymbols:>sign]

Coherent states with (Subscript[G, 1]^((2))) coupling, β=0.76 (left) and β=0.7 (right) with N=512

Integro-Differential Equation

Traveling wave ansatz

z(x,t)=(e^(-i Ω t))u(ξ≡x-ct)

The self consistency equation can be reduced into

c (Subscript[u, ξ])+i Ω u+(1/2)[(e^(-i α))U-(u^2)(e^(i α))(SuperStar[U])]=0,

where U≡K[u]. To solve the integro-differential equation:

1. Approximate u with Fourier modes

{(e^(-i n x)), (e^(-i(n-1)x)), ..., (e^(-i x)), 1, (e^(i x)), ..., (e^(i n x))}

# of real equations: 4n+2
# of real variables: 4n+2+2=4n+4
# of symmetries: 1+1=2 (translation and phase rotation)

2. Solve the truncated system by root finding algorithm and continue the system to get the solution branch {u(x,α),c(α),Ω(α)}.

Implementation

Load the solution obtained from simulation

ClearAll["Global`*"];
Get["initial.txt"];
Plot[tsol[x],{x,0,2Pi},AxesLabel->{Style["x",20],Style["θ",20]},PlotStyle->Black]

Approximate the solution with Fourier modes (n=64)

n=64;
z=Exp[I tsol[#]]&;
{a,b}=NIntegrate[z[x]E^(-I # x),{x,-Pi,Pi}]&/@{-1,1};
{ϕ,s}=-{#[[1]]+#[[2]],#[[1]]-#[[2]]}/2&@Arg[{a,b}];
fcoef=ParallelMap[Quiet[NIntegrate[z[x-s]E^(I(ϕ-# x)),{x,-Pi,Pi}]/(2Pi)]&,Range[-n,n]];
fcoef[[{n,n+2}]]=Re[fcoef[[{n,n+2}]]];

Project the integro-differential equation to truncated Fourier basis

n=64;
ai[1]=ai[-1]=0;
u=Total[(ar[#]+I ai[#])Exp[I # x]&/@Range[-n,n]];
U=(E^(-I x)) ar[-1]+(E^(I x)) ar[1]+(E^(-2 I x)) (ar[-2]+I ai[-2])+(E^(2 I x)) (ar[2]+I ai[2]);
Uc=U/.Complex[a_,b_]:>Complex[a,-b];

Construct the linear and nonlinear terms

uhat=(ar[#]+I ai[#])&/@Range[-n,n];
linear=((E^(-2I α)))Join[#,FourierCoefficient[U,x,#]&/@Range[-2,2],#]&[ConstantArray[0,n-2]]+(2(E^(-I α))I )(c1 Range[-n,n]+ Ω1)uhat;
linear=Module[{temp},temp=ComplexExpand[#];
Simplify/@{temp/.{I __-> 0},Total[Cases[temp,I __]/I]}]&/@linear;

u2hat=With[{temp=Collect[Expand[u^2],E^Times[__,x]]/.{E^Times[a__,x]/;Abs[a]>(n+2)-> 0}},
Insert[Sort[Cases[temp,E^Times[a__,x]b__:> {a/I,b},Infinity],#1[[1]]<#2[[1]]&][[All,2]],temp/.{E^Times[__,x]-> 0},n+3]];

Build the full system of equations.

eqs=Thread[Flatten[linear-ParallelMap[Module[{temp},temp=ComplexExpand[#];
{temp/.{I __-> 0},Total[Cases[temp,I __]/I]}]&,Total[(u2hat[[(#+6);;#]]&/@{-1,-2,-4,-5})(FourierCoefficient[Uc,x,#]&/@{-2,-1,1,2})]]]==0];

eqs[[1]]

Solve the system with FindRoot and continue the solution with respect to α

data=Reap[Sow[{β,sol=FindRoot[eqs/.{α->(Pi/2-β)},Join[{{c1,c/Pi},{Ω1,Ω/Pi}},Thread[{ar/@Range[-n,n],Re@fcoef}],DeleteCases[Thread[{ai/@Range[-n,n],Im@fcoef}],{0,0}]]]}];
Do[Sow[{β,sol=FindRoot[eqs/.{α->(Pi/2-β)},List@@@sol]}];,{β, 0.748,0.64,-0.002}];
][[2,1]];

data1=Reap[Sow[{β,sol=FindRoot[eqs/.{α->(Pi/2-β)},Join[{{c1,c/Pi},{Ω1,Ω/Pi}},Thread[{ar/@Range[-n,n],Re@fcoef}],DeleteCases[Thread[{ai/@Range[-n,n],Im@fcoef}],{0,0}]]]}];
Do[Sow[{β,sol=FindRoot[eqs/.{α->(Pi/2-β)},List@@@sol]}];,{β, 0.752,0.764,0.002}];
][[2,1]];

Results

Dependence of Ω and c on β

Get["data.txt"];
GraphicsRow[{ListLinePlot[Thread[{data[[All,1]],2Pi data[[All,2,2,2]]}],AxesLabel->{Style["β",15],Style["Ω",15]},PlotStyle->Black],ListLinePlot[Thread[{data[[All,1]],2Pi data[[All,2,1,2]]}],AxesLabel->{Style["β",15],Style["c",15]},PlotStyle->Black]},ImageSize->700]

Solution profile at different values of β

Manipulate[Plot[Evaluate[Mod[Arg[u/.data[[ind,2]]]-1.1,2Pi]-Pi],{x,-Pi,Pi},AxesLabel->{Style["x",15],Style["θ",15]},PlotStyle->Black,PlotRange->{{-Pi,Pi},{-Pi,Pi}},Epilog->Text[Style["β = "<>ToString@SetPrecision[data[[ind,1]],3],15],{2.25,2.9}]],{{ind,1," "},1,Length[data],1},Initialization:>(
n=64;
ai[1]=ai[-1]=0;
u=Total[(ar[#]+I ai[#])Exp[I # x]&/@Range[-n,n]];
Get["data.txt"];),TrackedSymbols:>ind]

Final Remark

We studied the spatial profile and stability of nonlocal phase-coupled oscillator systems via continuation by using Mathematica.

The following states are identified:

Splay states

Multi-cluster chimera states

Traveling waves

Reference

Multicluster and traveling chimera states in nonlocal phase-coupled oscillators, Phys. Rev. E, 90, 022919, (2014)

Initialization

$HistoryLength=0;
Share[];
SetDirectory[NotebookDirectory[]];
<<"data.txt";
PhaseOscillators1=Function[{phinitial,beta,tspan,mode,dt},
Module[
{n=Length@phinitial,dx,xx,totalsteps,phistore,phitemp,cossin,eiaN,btemp,atemp,k1,k2,k3,k4,count=2,count1=0,m,rhs,rk4coef},

If[FractionalPart[totalsteps=(tspan[[-1]]-tspan[[1]])/dt]!=0,Print[];Abort[]];
m=Round[(tspan[[2]]-tspan[[1]])/dt];
dx=2.0Pi/n;
xx=N[Range[n]]dx;

phistore=ConstantArray[0.,{Length@tspan,n}];
phistore[[1]]=phitemp=phinitial;

cossin={Cos[#],Sin[#]}&[mode xx];
eiaN=-Exp[-I beta]I dx dt;

rk4coef=Developer`ToPackedArray[{1.,2.,2.,1.}/6];
rhs=(btemp=Exp[-I#];Im[((eiaN Dot[#,btemp]&/@cossin).cossin)Conjugate[btemp]])&;

Nest[(phitemp+=rk4coef.FoldList[rhs[phitemp+#1 #2] &,rhs[phitemp],{0.5,0.5,1.}];
If[Mod[++count1,m]== 0,phistore[[count++]]=phitemp;])&,Null,Round@totalsteps];

phistore=Mod[phistore,2Pi]
]];

PhaseOscillators2=Function[{phinitial,beta,tspan,mode,dt},
Module[
{n=Length@phinitial,dx,xx,totalsteps,phistore,phitemp,cossin,eiaN,btemp,atemp,k1,k2,k3,k4,count=2,count1=0,m,rhs,rk4coef},

If[FractionalPart[totalsteps=(tspan[[-1]]-tspan[[1]])/dt]!=0,Print[];Abort[]];
m=Round[(tspan[[2]]-tspan[[1]])/dt];
dx=2.0Pi/n;
xx=N[Range[n]]dx;

phistore=ConstantArray[0.,{Length@tspan,n}];
phistore[[1]]=phitemp=phinitial;

cossin={Cos[mode #],Sin[mode #],Cos[(mode+1)#],Sin[(mode+1)#]}&[xx];
eiaN=-Exp[-I beta]I dx dt;

rk4coef=Developer`ToPackedArray[{1.,2.,2.,1.}/6];
rhs=(btemp=Exp[-I#];Im[((eiaN Dot[#,btemp]&/@cossin).cossin)Conjugate[btemp]])&;

Nest[(phitemp+=rk4coef.FoldList[rhs[phitemp+#1 #2] &,rhs[phitemp],{0.5,0.5,1.}];
If[Mod[++count1,m]== 0,phistore[[count++]]=phitemp;])&,Null,Round@totalsteps];

phistore=Mod[phistore,2Pi]
]];
xx=2.Pi Range[512]/512;